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Abstract. The ideal MHD equations are a central model in astrophysics, and their solution relies upon 
stable numerical schemes. We present an implementation of a new method, which possesses excellent stability 
properties. Numerical tests demonstrate that the theoretical stability properties are valid in practice with 
negligible compromises to accuracy. The result is a highly robust scheme with state-of-the-art efficiency. 
The scheme's robustness is due to entropy stability, positivity and properly discretised Powell terms. The 
implementation takes the form of a modification of the MHD module in the FLASH code, an adaptive mesh 
refinement code. We compare the new scheme with the standard FLASH implementation for MHD. Results 
show comparable accuracy to standard FLASH with the Roe solver, but highly improved efficiency and 
stability, particularly for high Mach number ffows and low plasma /3. The tests include ID shock tubes, 
2D instabilities and highly supersonic, 3D turbulence. We consider turbulent ffows with RMS sonic Mach 
numbers up to 10, typical of gas ffows in the interstellar medium. We investigate both strong initial magnetic 
ffelds and magnetic ffeld ampliff cation by the turbulent dynamo from extremely high plasma /3. The energy 
spectra show a reasonable decrease in dissipation with grid reffnement, and at a resolution of 512^ grid cells 
we identify a narrow inertial range with the expected power-law scaling. The turbulent dynamo exhibits 
exponential growth of magnetic pressure, with the growth rate twice as high from solenoidal forcing than 
from compressive forcing. Two versions of the new scheme are presented, using relaxation-based 3-wave and 
5-wave approximate Riemann solvers, respectively. The 5-wave solver is more accurate in some cases, and 
its computational cost is close to the 3-wave solver. 



1. Introduction 

In order to model complex nonlinear astrophysical phenomena, numerical solutions of the equations of 
ideal MHD are central. Examples include the dynamics of stellar atmospheres, star formation and accretion 
discs. Magnetised gas flows in astrophysics are typically highly compressible, nonlinear and often supersonic. 
Hence, obtaining stable numerical results in such regimes is extremely challenging if the numerical scheme 
must also be both accurate and efficient. Here we present a new MHD solver that preserves stable, physical 
solutions to the compressible MHD equations by construction, while at the same time showing improved 
efficiency and comparable accuracy to standard schemes based on the very accurate, but less stable Roe 
approximate Riemann solver [59 . 
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We consider numerical solutions of the ideal MHD system (here written in conservation form, and in three 
dimensions, letting Is denote the 3x3 identity matrix) 

+ V • (pu) = 0, 
{pu)t + V • (pu u + (p + ^ |Bp) /3 - B B) = 0, 
(1.1) + V • [(^ +p + ^|Bp)u - (B • u)B] = 0, 

Bt + V • (B (g) u - u (g) B) = 0, 
V-B = 0, 

with the specific internal energy e, such that the total energy density is given hy E = pe-\- ^pu^ + ^|Bp, and 
the pressure given by the equation of state p = p{p^ e). The Cartesian components of the velocity are denoted 
by u = {u^v^w). The system fits the generic form of a conservation law [/^ + V • F(/7) = 0, except for the 
restriction on V • B. However, if this restriction is satisfied at the initial time t = 0, it automatically holds 
at later times t > for the exact solution. Shock conditions are common in astrophysics, hence ( |1.1[ ) should 
be understood in a weak sense. The lack of regularity and large scale ranges in astrophysical flows make 
it particularly challenging to devise numerically stable schemes. Recent developments in nonlinear stability 
analysis have made such schemes possible, however. The central stability notions are entropy stability and 
posit ivity of mass density p and internal energy pe. This paper presents an implementation of the robust 
positive second order scheme from [72] , which uses the entropy stable approximate Riemann solvers of [H [5] . 

Several codes are available for high-performance compressible MHD simulations. We mention Athena [6F, 
AstroBEAR [13 , Chombo [60,, Enzo [73 , FLASH [23,, Nirvana [75 , Pluto [49 , RAMSES [22 and VAC 
[68] . The positive scheme of [72] has so far been unavailable in such codes. We have implemented it as an 
alternative MHD module in the FLASH code. 

We focus here on a comparison with the original FLASH code. Benchmark tests in one and two spatial 
dimensions are presented. Transition to unstable turbulence-like flow is particularly emphasised. Additional 
tests of the algorithms may be found in [5l [35] and [72^ . As an example application we consider turbulence in 
molecular clouds, which is characterised by high Mach numbers and the presence of magnetic flelds. Before 
presenting the results, we give a review of the numerical schemes and their theoretical properties. 

2. Numerical schemes 

In this section we give a quick summary of the numerical techniques used, and their justiflcation. Consider 
flrst a system of conservation laws in one spatial dimension (we consider the x-dimension) Ut -\- F{U)x = 0. 
The numerical schemes we will use are of the flnite volume type 

(2.1) 5Atf/, = C/,-^(7-,+ i 

where Ui are averages over intervals (or 'cells') of length h indexed by z at a given time t, and the operator 
Sai updates the cell averages to time t + At. The numerical fluxes T are evaluated at the cell interfaces, 
hence U is conserved and we call the scheme conservative. 

2.1. First order accuracy in one space dimension. First order accurate schemes can be given by 

(2.2) =^(C/i,C/i+i). 

Typically, is given by an approximate Riemann solver [3l[671|41]. The following are considered here: 

The Roe solver and the HLLE (Harten, Lax, van Leer and Einfeldt) solver as implemented in FLASH, 
and the multi-wave HLL-type solvers of [H [5] (named HLL3R and HLL5R, or HLLxR). Other HLL-type 
solvers are given in j27l |50l |26l Sll SH [25] among others. The numerical flux should be consistent (i.e. 
J-'{U,U) = F ([/)), and satisfy appropriate stability criteria. For gas and plasma dynamical systems a strong 
and physically meaningful stability criterion is implied by the second law of thermodynamics, which reads 
in one dimension: 

(2.3) {ps)t + {pus), > 0, 

s being the speciflc entropy of the system. This inequality holds, by deflnition, for an exact Riemann solver, 
and should also hold for an approximate Riemann solver. The HLLxR solvers were proved in [4, 5 to satisfy 



(2.3), and consequently the resulting conservative scheme satisfles a discrete version of (2.3). We refer to an 



approximate Riemann solver satisfying (2.3) as being entropy stable. For the Roe and HLLE solvers there 



are no rigorous proofs of entropy stability. For HLLE this appears not to be a problem, while Roe solvers 
yield non-physical shocks at sonic points. Therefore, Roe solvers are typically equipped with a so-called 
entropy-fix, which adds extra dissipation at sonic points. The Roe solver also tends to produce negative 
mass density values, while HLLE and HLLxR provably preserve positive mass density. Positivity of mass 
density and entropy stability imply that the internal energy remain positive. The positivity of density and 
internal energy ensures that the hyperbolicity of the equations does not break down (in which case the 
numerical schemes no longer make sense), and implies that the computed conserved quantities are stable in 
[55]. Entropy stability and positivity are referred to as nonlinear stability conditions since they apply 
directly to nonlinear systems. The 'R' in HLLxR comes from the derivation of the solvers from a relaxation 
approximation to the MHD equations, hence they are referred to as relaxation solvers. The underlying 
relaxation approximation has formally been found to be entropy dissipative through a Chapman-Enskog 
analysis [4 . The corresponding relaxation solver for hydrodynamics [2 , was tested in an astrophysical 
turbulence setting in [34j . 

Approximate Riemann solvers tend to vary in accuracy according to the level of detail they take into 
account. Typically, some wave modes are better resolved than others. For the solvers considered here, the 
HLLE solver can only optimally resolve the fast waves, while the Roe solver can, in principle, resolve all 
isolated waves as well as the exact solver. The relaxation solvers lie somewhere in between: The HLL3R can 
optimally resolve material contacts and tangential discontinuities at vanishing Alfven speed, and the HLL5R 
additionally resolves velocity shears at vanishing Alfven speed. Hence HLL5R is a true generalisation of the 
HLLC solver for the Euler equations. These differences in resolution most prominently appear when the 
wave modes in question are close to stationary, i.e. much slower than h/ At. 

2.2. Second order accuracy in one space dimension. In order to obtain second-order accuracy, some 
interpolation, i.e. reconstruction of the states from cell averages, has to be employed. We consider here 
the conservative MUSCL-Hancock scheme, introduced in ^ (see also [67 ). We refer to ^ for the non- 
conservative case. Let W denote the primitive state variables (p, u, B,p). Away from discontinuities the 
equations can be rewritten as Wt + A(W)iyaj = for a matrix ^4(1^). The MUSCL-Hancock algorithm goes 
as follows: 

(1) Reconstruction: evaluate discrete differences DWi. For oscillation control, we use the MC-limiter 
(monotonised central limiter), so for each component of Wi we take 

(2.4) DWi = cT^min (2\W,+, - W,\,^\W,+, - W,.,\,2\Wi - Wi-i 



with 

(2.5) (Ji = 

(2) Prediction step: evaluate 

(2.6) W,' = W,- ^A{W,)DWi 



1, W,^i -W,>0,W,- W,^i > 
-1, Wi^i -W^<0,Wi- Wi-i < 
0, otherwise. 



2h 



(3) Evaluate the cell edge values 



(2.7) Wr =W,'- \dW,, = W^^ ^DWi. 

(4) Use the cell edge values as input to the numerical fiux in the conservative scheme, 

(2.8) SAtUi = Ui-^{^ (C/+, - ^ (C/+ 1, U-)) , 

where U denotes the conserved state variables. 

Using primitive variables as the basis for reconstructing the states was recommended for example in [12]. 
It ensures that material contact discontinuities are reproduced exactly (and also shear waves at vanishing 
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Alfven speed). Alternatively, one can use that the primitive form Wt + = may be diagonalised 

as 



(2.9) {Ri)t^X^{Ri), = 0, j = l,2,...,d 

with Wi = XiRi. The matrix Xi is given by the eigenvectors X^ of A{Wi). The reconstructed gradients 
may then be evaluated using the relation DWi = XiDRi^ given gradients DRi^ and we get 

(2.10) wt = XiRf = Xi{Ri - ^XRi ± ^DRi) = Wi- ^A{Wi)DWi ± ^DWi, 

where A is the diagonal matrix having the eigenvalues A-^ as entries. The original FLASH code uses the 
characteristic variable-based reconstruction, while we use the primitive variables with the relaxation solvers. 



The limiter (2.4) forces the scheme to behave like the first order scheme near a discontinuity, hence 
ensuring that the entropy dissipation of the approximate Riemann solver kicks in. Even so, it turns out to 
be important to ensure the positivity of mass density and internal energy. In [72], a limiting procedure was 
presented such that the MUSCL-Hancock scheme is positive for any positive approximate Riemann solver. 
It consists of a modification that applies to any piecewise linear reconstruction of the W or R -variables. 
This limiter (more specifically the W-reconstruction' of [72]) is used with the relaxation solvers. 

2.3. Multidimensionality. Clearly our scheme is based on the treatment of one-dimensional Riemann 
problems. Multidimensional Riemann problems are much too complicated to be used as a basis of schemes, 
so various means have been suggested to derive multidimensional schemes from one-dimensional schemes. 
We only consider Cartesian grids here, where one can rely on dimensional splitting. This method consists 
of first solving the system ignoring the y and z- derivatives, then ignoring the x- and ^-derivatives etc. The 
splitting adds an error, and in order to maintain second order accuracy, the order of the directions must be 
reversed between each time step (Strang splitting). Unsplit methods have been developed to eliminate this 
error completely (e.g. in [39^, 'SS^). We limit our discussion to the split scheme, and focus on issues specific 
to the multidimensional MHD equations. Note however that these issues, as well as the proposed solutions, 
are not specific to split schemes. 

Due to considering the directions separately, we have to make sense of one-dimensional problems where the 
longitudinal magnetic field component is non-constant, in other words we can no longer assume V • B = 
initially. In order to handle the resulting more general one-dimensional problems, we use the following 
approach of [56 , consisting of modifying the evolution equation for B to 

^2 Bt + V • (B u - u B) - uV • B = 0. 



Equation (2.11) implies that 

(2.12) (V-B), + V-(uV-B) = 0, 

consequently V • B remains zero analytically, and numerical errors in V • B are transported by the flow. 
In the original FLASH code, the momentum and energy equations were also modified, which is what was 
actually recommended in [56 . We found it sufficient for stability in [5^ and [72 to only include the Powell 
term for the induction equation. 

The one-dimensional Powell system is only conservative when {Bn) is constant, so the schemes must take 
the form 

(2.13) S^tUi = U^-^ {HUt. U-^i) - MUt-i,U-)) - AtS,{U,-i,Ui, U,+i). 

In the original FLASH code, the source contributions from the cell edges are ignored, and the Powell source 
term is simply evaluated by a central discretisation. In [5 and [72 , we found that it was essential for 
stability, both practically and theoretically, to discretise the Powell term in a proper upwind manner. This 
was carried out in [4 for the first order scheme by extending the relaxation system to the Powell system, see 
also [^ . The second order extension can be found in [72] , and numerical tests demonstrating the importance 
of including a properly discretised Powell source are in [25", "SF and [72 . 

We include divergence cleaning in the form of the parabolic cleaning method of [47 ; see also [14 . We 
use it as already implemented in FLASH. It has the advantage of not expanding the numerical stencil of the 
scheme. FLASH also has the non-local projection method of [6 implemented, which could in principle be used 



Scheme HLL3R/HLL5R FLASH-Roe FLASH-HLLE 

Riemann solver HLL3R/HLL5R Roe HLLE 

Reconstruction positive, primitive characteristic characteristic 

Powell term B only full full 

Powell term discretisation upwind central central 
Table 1. Summary of the schemes. 



also with the new schemes. An alternative to using (2.11) is the staggered mesh (or constrained transport) 



approach (reviewed in [69]), which in one dimension essentially means evaluating Bn at the cell interface 
instead of as a cell average. The numerical stability results we rely on (entropy stability and positivity) 
are not available for staggered methods, but staggered schemes have the advantage of guaranteeing that 
V • B = to approximation order in smooth regions. Note that all approximate Riemann solvers we consider 
may also be used as components of staggered methods. 

2.4. Isothermal gas. As a simple model of cooling (e.g., due to radiation), it is common to assume that 
the gas temperature is constant. This means that p = c^p, where the constant Cg is the sound speed. A 
practical way to implement isothermality consists of (i) solving the full equations for ideal gas with a 7 
close to 1 within each time step, and (ii) projecting to the isothermal state between each time step. We 
use 7 = 1.00001. This projection method is preferred due its simplicity, accuracy, and generalisability. By 
generalisability we mean that the projection top = c^p is easily replaced with more elaborate cooling models. 
The projection method requires robust underlying numerics, since the full energy equation has to be solved 
properly. 

2.5. The schemes. The changes made to the original FLASH code in order to implement the new scheme 
are as described above. The time step size is chosen from the maximum wave speed of the previous time step 

and the maximum of + ^ over all cell averages (cg denoting sound speed). We use a CFL-number of 0.8 
throughout, unless otherwise noted. We used version 2.5 of FLASH, but it is straightforward to implement 
the new scheme in other finite volume codes. 

For numerical testing, we compare the new schemes HLL3R and HLL5R with the basic FLASH schemes 
for MHD. FLASH offers both a Roe and an HLLE approximate Riemann solver. The relevant schemes are 
summarised in Table [l] All schemes are finite volume schemes of MUSCL-Hancock type, with differences in 
the (approximate) Riemann solver, reconstruction method and Powell terms as described above. The names 
of the schemes are taken from the Riemann solvers for simplicity, although other algorithmic differences may 
also be important. To distinguish the schemes from the Riemann solvers, we use the 'sans serif font to 
indicate the scheme and normal font type for the Riemann solvers (see Tab. [T]). The thermodynamics are 
either adiabatic with an ideal gas equation of state or isothermal. 

3. Numerical studies ID 

In one spatial dimension, the differences between the schemes are given mainly by the accuracy of the 
approximate Riemann solver and the stability properties of both the approximate Riemann solver and the 
reconstruction. 

3.1. Brio-Wu shock tube. We first consider initial shock tube data from [8 that have become a standard 
test case. It contains a variety of numerically challenging wave types occurring in MHD fiows, in particular 
a fast rarefaction, a compound wave, a material contact, and a slow shock. The initial data are given by 
U = Ui for X < 0.5, and U = [/^ for x > 0.5, with 7 = 5/3 and 

pi = l,u^ =0,B^ = (0.75, 1,0), = 1 
Pr = 0.125, = 0, = (0.75, -1, 0),^^ = 0.1. 



Figure [33] shows small differences between HLL3R and FLASH-Roe on this standard test with FLASH-Roe 
slightly sharper on the compound wave. FLASH-HLLE is similar to HLL3R, except that the material contact 



was more smeared out than in HLL3R (see Fig. 3.2). We conclude that all the schemes perform reasonably 
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Figure 3.1. Brio-Wu shock tube. Resolution h = 256"^, time t = 0.1. FLASH-Roe and 
HLL3R are compared. The waves seen are from left to right, a fast rarefaction, a slow mode 
compound wave, a material contact discontinuity and a slow shock. The smaller amplitude 
right-going Alfven and fast waves are omitted for clarity. 




Figure 3.2. The same experiment as Figure 3.1 FLASH-HLLE and HLL3R are compared. 
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Figure 3.3. Low (3 expansion wave. Resolution h = 256 ^. Adiabatic case, velocity u = 3.1 
at time t = 0.05. FLASH-Roe and HLL3R are compared. 




Figure 3.4. Low f3 expansion wave, isothermal gas. Resolution h = 256 ^. Velocity 
u = 3.1 at time t = 0.05. FLASH-HLLE and HLL3R are compared. 



well on this test case, with some differences in accuracy depending on the level of detail in the approximate 
Riemann solver. 

3.2. Low 13 expansion waves. This is a test from [5 . It consists of rarefactions into a region of low plasma 
f3 (defined as /3 = ^). The initial data are given by U = Ui for x < 0.5, and U = Ur ioi x > 0.5, with 
7 = 5/3 and 

pi = l,uz = (-i/,0,0),B^ = (l,0.5,0),p^ = 0.45 
Pr = 1, = {u, 0, 0), = (1, 0.5, 0),pr = 0.45. 

Hence the sound speed is initially a/3/2 ^ 0.87. At = 3.1 there are only small differences between 
FLASH-Roe and HLL3R (Fig. [33]). 

We also performed this test for the isothermal case with = 0.45. FLASH-Roe could only handle u- 
values up to about 2.9. For higher values it failed to run due to negative pressure and density. In contrast. 



FLASH-HLLE and HLL3R proved quite stable also in the isothermal case. They are compared in Figure 3.4 
HLL3R is able to represent lower densities than FLASH-HLLE during the expansion. 



4. Numerical studies 2D 

In two dimensions we consider two standard test cases from the numerical literature and a Kelvin- 
Helmholtz instability. The tests exhibit strong shock waves and complicated wave interactions. Regarding 
numerical stability, in addition to the approximate Riemann solver and the reconstruction, the infiuence of 
different discretisations of the Powell term are investigated. 
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Figure 4.1. Orszag-Tang test. Pressure at time t = 0.5, with resolution h = 256 ^. Left: 
FLASH-HLLE, right: HLL3R. The result from FLASH-Roe was visually hard to discern from 
HLL3R. 




Figure 4.2. along a slice in the ^/-direction for the Orszag-Tang test. HLL3R (dashed) 
and FLASH-HLLE (solid). The result from FLASH-Roe was indistinguishable from HLL3R. 



4.1. Orszag-Tang. This is a standard test case for ideal MHD schemes. It consists of smooth velocity 
and magnetic field profiles going through shock formations, complicated wave interactions and eventually 
forming instabilities. The initial data are given by 7 = 5/3, 



(4.1) 



(p, u, B,p) = (1, sin(7r?/), sin(7rx), 0, — sin(7r?/), sin(27rx), 0, 0.6) . 



At t = 0.5 we observed only minor differences between HLL3R and FLASH-Roe, while FLASH-HLLE gave a 
smoother pressure peak along the central current sheet (see Fig. 4.1 ). Plotting Bx along the cut at x = 0.428 
(Fig. 4.2) shows that FLASH-HLLE gives a slightly more smeared-out current sheet. HLL3R and FLASH-Roe 
gave very similar results att = 0.5. Att = la magnetic island can be observed at the domain centre 
with FLASH-Roe and HLL3R, while with FLASH-HLLE this instability is suppressed by numerical viscosity 
(Fig. |4.3[ ). The HLL5R was also tested and gave results very similar to HLL3R and FLASH-Roe. 

We consider computational times in Table [2] They are from a computation on two processors, in the form 
of total CPU time in seconds and in the percentage of that computation time spent in the sweep routine 
(which is the only differing one). This should be taken with a grain of salt, as these numbers vary a bit 
from time to time. We took the average of two runs. It appears that HLL3R gives about a 10% speed-up 
compared to FLASH-Roe. FLASH-HLLE is also a little slower than HLL3R, which is due to the different 
reconstruction procedures. All in all we would expect FLASH-HLLE to be faster than FLASH-Roe, since only 
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Scheme HLL3R HLL5R FLASH-HLLE FLASH-Roe 

N time steps 477 539 476 482 

Sweep time fraction 43.6% 47.6% 47.2% 50.0% 
Average time / step 0.555 s 0.520 s 0.571s 0.604 s 

Table 2. Computational times for the Orszag-Tang test. 



the computation of the numerical flux differs. The time to compute the numerical flux in HLL3R should 
be somewhere in between the time for FLASH-Roe and FLASH-HLLE, and in addition HLL3R has a simpler 
reconstruction and a slightly more elaborate Powell term than FLASH-Roe and FLASH-HLLE. The HLL5R 
scheme is the most efficient, however it requires slightly smaller time steps. The shorter time steps are not 
surprising, since the wave speeds of the HLL5R Riemann solver may be up to 10% larger than with the 
corresponding 3- wave solver in some configurations; see ^. Consequently, there would be a roughly 10% 
reduction in time step size. 

4.2. Cloud-shock interaction. The initial conditions are the same as in [25 . The computational domain 
is (x, y) G [0, 2] X [0, 1] with artificial Neumann type boundary conditions. The initial conditions consist of 
a shock moving to the right, located at x = 0.05, and a circular cloud of density p = 10 and radius r = 0.15 
centred at (x^y) = (0.25,0.5). We took 7 = 5/3 and 

'3.86859 ifx< 0.05, 
(4.2) p = { 10.0 if {x - 0.25)2 + (y - 0.5)^ < 0.15^ 

1.0 otherwise, 



(4.3) 

(4.4) p = 

(4.5) B = 



(11.2536,0,0) ifx < 0.05, 
(1.0,0,0) otherwise, 

'167.345 if X < 0.05, 
1.0 otherwise, 

'(0, 2.18261820, -2.18261820) if x < 0.05, 
(0, 0.56418958, 0.56418958) otherwise. 



Hence, the shock has a sonic Mach number around 11.8 and an Alfvenic Mach number around 19.0. The 
solution becomes very complicated as the shock interacts with the cloud, and instabilities are observed 
along filaments of the cloud. We used these initial data to test the ability of the new scheme to handle 
adaptive mesh refinement. Eight levels of refinement were used with the highest resolution corresponding to 
a resolution of 1024 x 2048 grid cells. The refinement criterion was based on the mass density p with the 
FLASH code parameters ref ine_cutof f = 0.75 and deref ine_cutof f = 0.25. Mass density profiles from 



HLL3R and FLASH-Roe are shown in Figure 4.4 The front of the cloud looks more rounded with HLL3R, and 
the structures of the turbulent wake differ somewhat. Both FLASH-Roe and HLL3R agree qualitatively well, 
but since there is no analytic solution it is unclear how to make a quantitative assessment of the differences. 
HLL5R gave very similar results to HLL3R in this test. FLASH-HLLE roughly reproduced the shape of the 
cloud front, but gave a different structure in the turbulent wake. 

4.3. Kelvin-Helmholtz instability. The Kelvin-Helmholtz instability is important as a model for the 
onset of turbulence. From a numerical point of view, we found it illustrative of some differences between 
the Riemann solvers. We consider a case where both the velocity and magnetic fields switch polarity along 
a straight line. The initial data are given hy p = 3/5, 7 = 5/3, {v, By^ Bz) = and 

(4.6) p = 2,u = 0.5,5^ 0.5, 
for ^ > 0, and 

(4.7) p=l,u = -0.5, B^ = -0.5, 
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Figure 4.4. Cloud-shock interaction with 8 levels of refinement, corresponding to 1024 x 
2048 cells. Grid and mass density. Top: FLASH-Roe, bottom: HLL3R 

for y < 0. We perturb the velocity in the ^/-direction by 

(4.8) V = 0.0125 sin(27rx) exp [-I00{y - 0.5)^] . 

At y = —0.5 and y = 0.5 we impose reflecting boundary conditions. Hence, we consider a small perturbation 
to a tangential discontinuity in both u and B along the line ^ = 0. Initially the relative sonic Mach number 
is 1, the magnetosonic Mach number is 2/>/5, while plasma /3 = 87. We used a resolution of 512^ cells. 
In numerical experiments of this type, instability growth tends to occur at the small scales, and therefore 
depend on numerical viscosity. In [34 , this was reflected in a strong dependence between grid resolution and 
instability growth rate. Mass density profiles at time t = 1 are shown in Figure |4.5| The higher numerical 
viscosity of the HLLE solver results in FLASH-HLLE suppressing all instabilities at this time. HLL3R yields 
a sharper shear-discontinuity than FLASH-HLLE, but the roll-ups still did not develop. In contrast, HLL5R 
and FLASH-Roe produce well-developed Kelvin-Helmholtz roll-ups. 

The different behaviour of the schemes can be explained by the approximate Riemann solvers. By the 
remarks at the end of section |2.1[ Roe and IILL5R exactly resolve the unperturbed initial configuration, 
while IILL3R smears out the velocity shear slightly. HLLE smears out the material and magnetic contact 
discontinuities of the initial configuration, with the result that the instability is entirely suppressed at early 
times. The presence of small perturbations to grid-aligned stationary discontinuities makes this test case 
particularly sensitive to these properties of the Riemann solvers. In particular, we see the effect of exactly 
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Figure 4.5. Kelvin-Helmholtz instability at time t = 1. Top to bottom: FLASH-Roe, 
HLL5R, HLL3R, FLASH-HLLE. 
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Figure 4.6. Kelvin-Helmholtz at time t = 9, FLASH-Roe (left) and HLL5R (right) 



resolving the discontinuities, i.e., having a scheme that is 'well-balanced' with respect to the unperturbed 
state. In non-stationary or non-grid-aligned cases, none of the schemes have this well-balancing property, 
hence we expect differences to be less pronounced. 

At late times of the instability development, it is harder to assess the performance of the schemes. Results 
from FLASH-Roe and HLL5R are shown in Figure 4.6 Both schemes give a similar structure of two larger 



vortices. FLASH-Roe produces small-scale waves in the surrounding medium, which are possibly a numerical 
artefact. One possible source of spurious features in FLASH-Roe is the central discretisation of the Powell 
source term. Such errors have been pointed out in [25] , [24] , [35 and [72 . 

Finally, in order to demonstrate the robustness of HLL3R and HLL5R, we tested an isothermal case with 
a higher sonic Mach number. The sound speed was set to 0.1, and the initial data were {v^w^By^Bz) = 
everywhere, 

0.5,5, 



1, u 



2, 



for ^ > 0, and 



-0.5,5, 



for ^ < 0. The perturbation was again given by (4.8). Hence, the relative sonic Mach number was 10, and 



the relative Alfvenic Mach number was 1/a/2. Resolution was uniform with 512^ grid cells. The standard 
schemes, FLASH-Roe and FLASH-HLLE failed on this case due to unphysical states, but the new schemes. 



HLL3R and HLL5R produced physical results. Density profiles from HLL3R and HLL5R are shown in Fig. 4.7 



The resulting instability contains vortex-like structures with complicated internal features. Due to the chaotic 
nature of these solutions, we cannot evaluate the differences seen between the schemes. 
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Figure 4.7. High sonic Mach number Kelvin-Helmholtz test for isothermal gas at time 
t = 1. Top: HLL3R. Bottom: HLL5R. 
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5. Supersonic MHD turbulence 



In this section we investigate the performance of three of the MHD schemes, FLASH-Roe, HLL3R and 
HLL5R on the modehing of supersonic MHD turbulence. Testing the schemes in the regime of highly com- 
pressible turbulent flows is of particular interest for astrophysical applications ranging from molecular cloud 
formation on galactic scales down to scales of individual star-forming accretion disks. These disks can launch 
powerful supersonic jets, driven by magneto-centrifugal forces [1 and magnetic towers [45 . Here however, 
we concentrate on intermediate scales of molecular cloud evolution, which is largely determined by the inter- 
play of supersonic turbulence and magnetic fields. Understanding the nature and origin of supersonic MHD 
turbulence is considered to be the key to understanding star formation in molecular clouds [46l [151 EB SH] • 



5.1. Numerical setup and post-processing. We start with a short description of the type of setup often 
adopted in numerical models of molecular cloud turbulence. It is called turbulence-in-a-box, because a three- 
dimensional box with periodic boundary conditions and random forcing is used to excite turbulence. We use 
an isothermal equation of state throughout all turbulence-in-a-box runs, i.e., the pressure, pth, is related to the 
density by pth = with the constant sound speed, Cg. For these tests we used an initially uniform density, 
(p) = 1.93 X 10"^-^ gcm~^, corresponding to a mean particle number density of (n) = (p)/(/imp) ^ 500 cm~^, 
with the mean molecular weight, /i = 2.3 and the proton mass, mp, and a box size of = (4pc)^. The 
sound speed for the typical gas temperatures in molecular clouds (about 11 K) is roughly Cg = 0.2kms~^. 
For driving the turbulence we use the same forcing module as used in [TTl [181 EH ISH] , which is based on 
the stochastic Ornstein-Uhlenbeck process [16 . The spectrum of the forcing amplitude is a paraboloid in 
Fourier space in a small wavenumber range 1 < |k| < 3, peaking at |k| = 2. This corresponds to injecting 
most of the kinetic energy on scales of half of the box size. The turbulence forcing procedure is explained 
in detail in [19 . We use a solenoidal (divergence-free) forcing throughout this study unless noted otherwise. 
For the comparison of the different MHD schemes the forcing amplitude was set such that the turbulence 
reaches a root-mean-squared (RMS) Mach number, ^ 2 in the fully developed, statistical steady regime. 
The resolution was fixed at 256^ grid cells. We also consider a highly supersonic case with ~ 10 in 
section [5^ for which we investigate the resolution dependence by comparing runs with 128^, 256^ and 512^ 
grid zones. This study, however, could only be followed with the new types of solvers presented here (HLL3R 
and HLL5R), because the Roe solver is numerically unstable for highly supersonic turbulence. As shown in 
previous studies [TS] |64l |19] turbulence is fully developed after about two large-scale eddy turnover times, 
2T, where one large-scale eddy turnover time is defined as T = L/{2M.c^). The RMS velocity is thus defined 
as y = L/(2T) =Mc,. 

In order to test the MHD schemes in a strongly magnetised case, we set the initial uniform magnetic field 
strength to Bq = 8.8 x 10^ Gauss in the z-direction, which gives an initial plasma beta, /3 = pth/Pm,o ^ 0.25, 
with the magnetic pressure, Pm, o = ^^/{STr). No initial fiuctuations of the magnetic field were added, 
such that the initial RMS magnetic field is equal to Bq. However, the turbulence twists, stretches, and 
folds the initially uniform magnetic field, such that the RMS field strength increases due to turbulent 
dynam o am plification until saturation (see the review by Brandenburg and Subramanian 2005 [7 ; and 
section 



5.7). The Alfvenic Mach number is A4a = V/va = A4^y/3/2 ^ 0.7, i.e., we consider slightly sub- 



Alfvenic turbulence. The HLL3R scheme, however, was used in both the highly sub- Alfvenic {A4a <^ 1) 
and the highly super- Alfvenic (A^a ^ 1) turbulent regimes by [iniE], who modelled turbulent flows with 
f3 = 0.01 ... 10, Ma = 0.4 ... 50 and M = 2 .. . 20. 

All results were averaged over one eddy turnover time, T, in the regime of fully developed turbulence, 
t > 2T [19 . Since we produced output files every 0.1 T, we computed Fourier spectra and probability 
distribution functions for each individual snapshot and averaged afterwards over 11 snapshots for 2.3 < 
t/T < 3.3. The averaging procedure allows us to study the temporal fiuctuations of all statistical measures 
and to improve the statistical significance of the similarities and the differences seen between the different 
schemes. This is an important advantage over comparing instantaneous snapshots, because it must be 
kept in mind that individual snapshots can be different due to intermittent turbulent fiuctuations. These 
fiuctuations can either make the results look very similar or very different between the different schemes, 
if one compares instantaneous snapshots only. Thus, a meaningful comparison can only be made by time 
averaging the results and considering the statistically average behaviour of the schemes (see also, [58]). 
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Figure 5.1. Time evolution of the RMS Mach number, M for FLASH-Roe, HLL3R and 
HLL5R, shown as sohd, dotted and dashed hues, respectively. Turbulence is fully developed 
for t > 2T, where T is the large-scale eddy turnover time. 



5.2. Time evolution. Figure 5.1 shows the time evolution of the RMS Mach number, A4 = (u^)/cs. 
The magnetised fluid is accelerated from rest by the random stochastic forcing as discussed in the previous 
section. It reaches a statistical steady state at around t ^ 2T, such that turbulence is fully developed for 
t>2T (see also, [18l[19]). We can thus safely use time- averaging for 2.3 < t/T < 3.3 of the Fourier spectra 
discussed in the next section. 

5.3. Fourier spectra. As a typical measure used in turbulence analyses (e.g., [2TJ |34l [32] ) , we show Fourier 
spectra of the velocity and magnetic fields in Figure 5.2, decomposed into their rotational part, and 



their longitudinal part, EK For the velocity spectra (top panel), E^ ex |u^p is a measure for the rotational 
motions of the turbulent fiow, while oc |u|| p is a measure for the compressions induced by shocks. Since 
here we drive the turbulence with a solenoidal random force [19], and since the RMS Mach number is only 
slightly supersonic {A4 ~ 2) the bulk of the turbulent motions is in the rotational part, E^. From the 
velocity spectra we conclude that FLASH-Roe is the least dissipative scheme of the three schemes tested. 
The HLL3R is slightly more dissipative, which is seen from the faster loss of kinetic energy at wavenumbers 
k > 20. The HLL5R is in between the FLASH-Roe and the HLL3R spectrum. Given that the limited numerical 
resolution does not allow for an accurate determination of the turbulence spectrum in the inertial range, 
the differences for /c > 20 are acceptable (see also section [^6] below) . [32] and [19] find that about 30 grid 
cells are necessary to resolve the turbulence on the smallest scales. This means that the turbulence on all 
wavenumbers k > A^/30 ~ 8 may be affected by the limited resolution (here N = 256), which means that 
the differences between the different schemes for A: > 20 can be safely ignored. 



The magnetic field spectra (Figure 5.2, bottom panel) display small differences for A: > 10. As for 
the velocity spectra the differences occur on scales smaller than the scales that are affected by numerical 
resolution, which means that these differences are negligible in applications of supersonic MHD turbulence. 
The rotational part of the magnetic field, E^, shows slightly less dissipation for the HLL3R and HLL5R 
compared FLASH-Roe. The longitudinal part, £^m, is a measure for (V • B)^. The FLASH-Roe scheme 
keeps the divergence of the magnetic field smaller by an order of magnitude compared to HLL3R and the 
HLL5R on all scales. However, all schemes maintain the V • B constraint within acceptable values. The 
ratio J Em ^k/ J {E^ + E^) dk < 5.5 x 10~^ for all times, hence numerical V • B effects did not have any 
significant dynamical infiuence. 

5.4. Probability distribution functions. To further investigate the dissipative properties of the different 



schemes we show probability distribution functions (PDFs) of the vorticity, |V x u|, in Figure 5.3 (top left 
panel). The vorticity is a measure of rotation in the turbulent fiow. Stronger numerical dissipation leads to 
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Figure 5.2. Comparison of the Fourier spectra of the turbulent velocity {top panel) and 
the turbulent magnetic field {bottom panel) for the FLASH-Roe, HLL3R and HLL5R, shown 
as solid, dotted and dashed lines, respectively. The numerical resolution was 256^ grid cells. 
Both the velocity and the magnetic field spectra were decomposed into their rotational and 
compressible parts, and respectively. The temporal variations of the spectra are on 
the order of three times the line width in this plot. 



a faster decay of small-scale eddies. Figure |5.3| (top left panel) shows that the most dissipative among the 
three schemes tested, HLL3R produces smaller values of the mean vorticity by about 10% compared to the 
FLASH-Roe. The mean vorticity achieved with HLL5R is equivalent to FLASH-Roe. The RMS values of the 
vorticity for HLL3R and HLL5R are 11% and 1%, respectively smaller than for FLASH-Roe. This confirms 
the expected ranking of the schemes in terms of their numerical dissipation from the previous section: the 
HLL3R is the most dissipative, while the HLL5R is only marginally more dissipative than FLASH-Roe. 



The PDF of the divergence of the velocity field, V • u is shown in the top right panel of Figure 5.3 The 



asymmetry between rarefaction (V-u > 0) and compression (V-u < 0) is typical for compressible, supersonic 
turbulence (see, e.g., [54l|64]). FLASH-Roe produces slightly more zones with higher compression, which can 
be attributed to its lower numerical dissipation compared to HLL3R and HLL5R (cf., the Fourier spectra 
discussed in the previous subsection). However, all schemes agree to within the error bars, indicating strong 
temporal variations of compressed regions due to the intermittent nature of the shocks. 

The PDF of the gas density is an essential ingredient for analytic models of star formation (e.g., [52) l37t 
[29]). All those models need the density PDF to estimate the mass fraction above a given density threshold 
by integrating the PDF. Thus, many studies have focused on investigating the density PDF in a supersonic 
turbulent medium (e.g., I 111 |M1 IMl 133 El [13 Ell [40l [1^ Figure [53] (bottom left panel) shows the PDF 
of the logarithmic density, s = ln(p/(p)), where (p) is the mean volume density. In this transformation 
of the density, the PDF follows closely a Gaussian distribution, i.e., a log- normal distribution in p with 
a standard deviation of = In (l + PM'^) (see, e.g., ^9]). From this expression we derive a turbulence 
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Figure 5.3. Probability distribution functions of the vorticity, |V x u| {top left)^ the di- 
vergence of the velocity, V • u {top right)^ the logarithmic density, s = ln(p/(p)) {bottom 
left)^ and the magnetic pressure, Pm/Pm,o = ^^/^o {bottom right), for FLASH-Roe, HLL3R 
and HLL5R, shown as solid, dotted and dashed lines, respectively. The numerical resolution 
was 256^ grid cells. The grey hatched regions indicate the 1-a temporal fluctuations of the 
distributions. 

forcing parameter of 6 = 0.35 ±0.05 for the given RMS Mach number, M = 2.3 ±0.2 and ag = 0.71 ±0.04 for 
all three schemes. This is in very good agreement with the expected forcing parameter for purely solenoidal 
forcing as applied in this study, 6 = 1/3 (see [TT, 19 J. 

Finally, we show the PDF of the magnetic pressure, p^ = B'^ / {Sir) in units of the initial magnetic 
pressure, Pm,o, in Figure [53| (bottom right panel). No significant differences are seen between the three 
different schemes. Thus, the magnetic pressure distribution is quite robust. 

5.5. Computational performance of the MHD schemes. All runs with the three different schemes, 
FLASH-Roe, HLL3R and HLL5R were performed on the same machine, the HLRB II at the Leibniz Super- 
computing Center in Munich in a mode of parallel computation (MPI) with 64 CPUs. In Table [s] we list 
the wall clock time per step, the total number of steps, the total amount of CPU hours spent, and the CFL 
number used for each scheme in our driven MHD turbulence experiments with a numerical resolution of 
256^ grid cells. Only for FLASH-Roe we had to use a CFL number of 0.2 instead of 0.8 (both HLL3R and 
HLL5R), because it became unstable for CFL numbers > 0.2 and crashed due to the production of cells 
with negative densities. The time per step was estimated from averaging over 0.1 T within the last turnover 
time, i.e., between 3.2 < t/T < 3.3 by averaging over 309, 70 and 71 steps for the FLASH-Roe, HLL3R and 
HLL5R, respectively to avoid variations produced by I/O processing. Considering the time per step, HLL3R 
and HLL5R are about 4% and 3%, respectively faster than FLASH-Roe. This speed-up of HLL3R is slightly 
less than the speed-up measured for the Orszag-Tang test (cf. Tab.|2|. The difference between HLL3R and 
HLL5R is almost negligible. Indeed, the two corresponding solvers have a very similar structure, except that 
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Scheme 


Time per step [s] 


Total steps 


Total CPU hours 


CFL 


FLASH-Roe 


69.8 ±0.5 


16,502 


20,480 


0.2 


HLL3R 


67.3 ±0.5 


1,920 


2,297 


0.8 


HLL5R 


67.9 ±0.5 


2,065 


2,493 


0.8 



Table 3. Performance of the MHD schemes in driven MHD turbulence with = 2 for 
a numerical resolution of 256^ grid cells. All runs were performed in a mode of parallel 
computation with 64 CPUs on the SGI Altix 4700 (HLRB II) at the Leibniz Supercomputing 
Center in Munich. Note that the run with the FLASH-Roe scheme was unstable for CFL 
numbers > 0.2 in this test. 



the HLL5R solver typically goes through one more if-test compared to the HLL3R solver. HLL5R goes 
through a few more time steps, which is either because of its slightly more restrictive CFL constraint, or 
flow features such as local differences Alfven speed. 

The total number of steps and the total CPU time spent by tend = 3.3 T for each run are shown in the 
third and fourth column of Table [H Given the factor of 4 between the CFL number used with FLASH- 
Roe compared to HLL3R and HLL5R it is not surprising that the total number of steps is much higher for 
FLASH-Roe. However, the number of steps for FLASH-Roe is about 8.6 and 8.0 times higher compared to 
HLL3R and HLL5R respectively, which is more than the difference in the CFL numbers. The reason for this 
is that FLASH-Roe is close to being unstable even with CFL = 0.2. This produces timesteps that are often 
significantly smaller than necessary. Taken altogether, for the problem setup discussed in this section, HLL3R 
and HLL5R are factors of 8.9 and 8.2 respectively faster than FLASH-Roe, which is reflected by the total 
amount of CPU time spent for a successful MHD turbulence run with A4 ^ 2. For higher Mach number 
runs, however, FLASH-Roe was extremely unstable, which did not allow us to make this comparison of the 
schemes at higher Mach number. 



5.6. Resolution study. As outlined in the introduction, interstellar turbulence is highly supersonic with 
typical Mach numbers of the order of ~ 10. As discussed in the comparison of the different MHD 
schemes in the previous subsections, the FLASH-Roe scheme turned out to be highly unstable even for low 
Mach number turbulence. This limitation is overcome with the HLL3R and HLL5R schemes. In this section, 
we apply our new HLL3R scheme to MHD turbulence with = 10 as an example of the high Mach number 
turbulence typical of the interstellar medium and discuss numerical resolution issues. The three runs analysed 
in this section were run at numerical resolutions of 128^, 256^ and 512^ grid cells with a slightly smaller 
initial magnetic field, = 4.4 x 10~^ Gauss, which gives an initial plasma beta of /3 ~ 1. 

First, we show slices through the three-dimensional box at = L/2 of the density and magnetic pressure 



fields in Figure 5.4 With increasing resolution more small-scale structure is resolved. The density and 
magnetic pressure are correlated due to the compression of magnetic field lines in shocks. Vector fields 
showing the velocity and the magnetic field structure are overlaid on the slices of the density and magnetic 
pressure, respectively. 

Numerical dissipation becomes smaller with higher resolution. Since the inertial range of the turbulence 
is defined as all the length scales much smaller than the forcing scale and much larger than the dissipation 
scale, a minimum numerical resolution of about 512^ grid cells is required to separate the inertial range 
from the forcing and dissipative scales. To see how dissipation depends on numerical resolution we plot the 



velocity and magnetic field spectra in Figure |5.5| for the three different resolutions in analogy to Figure |5.2[ 
The turbulence is driven on large scales corresponding to k ^ 2 (injection scale, i^inj), while the spectra 
become steeper and begin dropping to zero at wavenumbers k > 10, 20, and 40 for 128^, 256^, and 512^, 
indicating the onset of dissipation (on scales ^dis)^ respectively. Thus, the range of scales separating the 
injection scale from the dissipation scale increases with increasing resolution. However, the inertial range 
is defined as all scales, i with Linj ^ i ^ idis- Purely hydrodynamical simulations of driven turbulence 
with up to 1024^ grid cells show that the inertial range just becomes apparent for resolutions > 512^, which 



means that we cannot see a convergence of the slope in the spectra shown in Figure 5.5 However, measuring 
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Figure 5.4. Slices of the density {left) and magnetic pressure {right) at y = L/2 and 
t = 2T for numerical resolutions of 128^ {top), 256^ {middle), and 512^ {bottom). 
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Figure 5.5. Comparison of the Fourier spectra of the turbulent velocity {top panel) and 
the turbulent magnetic field {bottom panel) for HLL3R in = 10, /3 = 1 MHD turbulence. 
The infiuence of the grid resolution is shown: 128^ (dotted), 256^ (dashed) and 512^ (sohd). 
Both the velocity and the magnetic field spectra were dec omp osed into their rotational and 
compressible parts, E-^ and respectively, as in Figure 5.2 



the power-law slope, E oc of the velocity spectra in the wavenumber range 5</c<10at512^ yields 
slopes of /3 « —1.6 ±0.1, and /3 ~ — 1.9±0.1, for E^ and ^Jl, respectively. The solenoidal part, E^ is 
in good agreement with the Kolmogorov spectrum {(3 = —5/3) [Ml EI], while the compressible part, E^ is 
significantly steeper and closer to Burgers turbulence {/3 = 2) [11], applicable to a shock-dominated medium. 
Note that E^ clearly decreases with increasing resolution. 

5.7. Turbulent dynamo. Magnetic fields are ubiquitous in molecular clouds, but it remains controversial 
whether these fields have an infiuence on the cloud dynamics (see, e.g., [46 ). However, it is widely accepted 
that magnetic fields play a significant role on the scales of protostellar cores, where they lead to the generation 
of spectacular jets and outflows, launched from the protostellar disks, a process for which a wound-up 
magnetic field seems to be the key [Ij [45]. Thus, we test in this section whether our new MHD schemes 
can represent wound-up magnetic field configurations with the same turbulence-in-a-box approach as in the 
previous subsections. 

The magnetic pressure can become comparable to the thermal pressure in dense cores due to the am- 
plification of the magnetic field through first, compression of magnetic field lines, and second, due to the 
winding, twisting and folding of the field lines by vorticity, a process called turbulent dynamo (see [7 for 
a comprehensive review of turbulent dynamo action in astrophysical systems). Magnetic field amplification 
in the early universe during the formation of the first stars and galaxies was discussed analytically in [63 . 
showing that the magnetic pressure can reach levels of about 20% of the thermal pressure in primordial 
mini-halos, thus potentially infiuencing the fragmentation of the gas. Numerical simulations also show that 
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Figure 5.6. Shows the magnetic pressure as a function of time (for 200 eddy turnover times, 
T) for solenoidal forcing (sohd hue) and compressive forcing (dashed hue). The dash-dotted 
horizontal hue shows the thermal pressure. The turbulent dynamo works with an exponential 
growth rate of about 0.60/T for solenoidal forcing and 0.28/T for compressive forcing (dotted 
lines). For comparison, typical amplification rates found in subsonic, solenoidally driven 
turbulence are about 0.5/T (e.g., [62]). 



the turbulent dynamo is an efficient process to amplify small seed magnetic fields during the formation of 
the first stars and galaxies (j74l[66]). In present-day star formation the infiuence of magnetic fields on the 
scales of dense cores were investigated in numerical work ( [57l ^0 ) , concluding that magnetic fields strongly 
affect fragmentation of dense gas. Understanding the magnetic field growth due to the turbulent dynamo 
is thus crucial for future studies of star formation. Many studies have focused on subsonic turbulence (e.g., 
[62 ), with only very few contributions on the supersonic regime. For instance, |28J only studied the turbulent 
dynamo for mildly supersonic Mach numbers, M. ^ 2.5. For molecular clouds however, the highly supersonic 
regime is more relevant. 

Figure 5.6 shows a comparison study of the turbulent dynamo operating in the supersonic regime {M. = 5). 
We used the HLL3R for this test with a numerical resolution of 128^ grid cells. Unlike the previous turbulence 
runs discussed above we started from an extremely small initial magnetic field, Bq = 4 Ax 10~^^ Gauss, which 
corresponds to an initial plasma beta of ^ ^ 10^^. We used two limiting cases to drive the turbulence: purely 
solenoidal (divergence-free) forcing as above, and purely compressive (curl- free) forcing (see [19] for a detailed 
analysis of the differences of solenoidal and compressive turbulence forcings and their role in the context 
of molecular cloud dynamics and star formation). Figure 5.6 shows that the turbulent dynamo leads to 
an exponential growth of the magnetic pressure over more than 15 orders of magnitude within about 100 
large-scale eddy turnover times, T. The dynamo growth rate is about twice as large for solenoidal forcing 
(^ 0.60/T) compared to compressive forcing 0.28/T), due to the higher average vorticity generated by 
solenoidal forcing (see [E]), which makes the small-scale dynamo more efficient. The actual growth rates, 
however, depend on the kinematic and magnetic Reynolds numbers. Since we did not add physical viscosity 
and magnetic diffusivity, these numbers are controlled by numerical viscosity and diffusivity, and thus by 
numerical resolution. A resolution study of the dynamo in recent FLASH simulations is presented in [66] 
and in [20 , the kinematic and magnetic Reynolds numbers found in dynamo simulations similar to the ones 
presented here are around 200 for a numerical resolution of 128^ grid cells. A more detailed analysis of the 
Mach number dependence of the turbulent dynamo amplification in solenoidal and compressive forcings is 
in preparation. The dynamo saturates after about TOT and 140 T for solenoidal and compressive forcing, 
respectively. The saturation level is close to the thermal pressure in both cases. For compressive forcing the 
saturated magnetic pressure is about 5% of the thermal pressure, while it is 40% for solenoidal forcing. This 
numerical test shows that the turbulent dynamo works with the new HLL3R scheme. This is an important 
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test as it shows that the scheme reproduces the expected amphfication of the magnetic pressure due to the 
winding-up of magnetic field fines. 

6. Summary 

We presented an implementation of an accurate, efficient and fiigfily stable numerical metfiod for MHD 
problems. Tfie metfiod is reviewed fiere, and presented in detail in It is implemented as a modification of 
tfie FLASH code [23], wfiicfi enables large-scale, multi-processor simulations and adaptive mesfi refinement. 
In [72] it was found tfiat our metfiod could fiandle significantly larger ranges of tfie sonic Macfi number and 
plasma /3 tfian a standard MHD scfieme. Tfiis was confirmed in tfiis paper by comparisons witfi tfie standard 
FLASH code. Tfie algoritfimic cfianges underlying tfie increased stability can be broken down into tfiree 
parts: 

(1) An entropy stable approximate Riemann solver tfiat preserves positivity of density and internal 
energy ([lll^). 

(2) For second order accuracy, a reconstruction metfiod tfiat ensures positivity (|72j). 

(3) In multidimensions, a stable discretisation of tfie Powell system ([77). 

All tfiese ingredients were essential in obtaining tfie desired stability and efficiency of tfie overall scfieme. Tfie 
different elements of tfie new scfieme fiave been studied separately in previous papers. Wfiile [72] focused 
on tfie positive second-order algoritfim and multidimensionality, only a single approximate Riemann solver, 
HLL3R, was considered. Tfie present study contrasts a standard scfieme to tfie combination of tfiese tfiree 
new ingredients. 

Tfie new scfieme was implemented in two versions featuring tfie 3-wave (HLL3R) and 5-wave (HLL5R) 
approximate Riemann solvers of [5 respectively, wfiile tfiere were two standard implementations of tfie 
FLASH code (version 2.5), using tfie Roe and tfie HLLE approximate Riemann solvers. We observed some 
increase of numerical dissipation compared to tfie Roe solver of FLASH, but it was minor, and due to tfie 
replacement of tfie Roe solver witfi tfie robust and efficient HLL-type solvers. Tfie HLLE solver was found 
to be tfie most dissipative, wfiile HLL5R sfiowed almost identical dissipation properties and accuracy to tfie 
Roe solver. HLL3R was ranked between HLLE and Roe in terms of accuracy. 

As a pfiysical application, we fiave considered forced MHD turbulence at fiigfi sonic Macfi number. We 
were able to compare tfie new and old scfiemes at RMS sonic Macfi number 2 witfi an initial plasma (3 = 0.25. 
Tfie scfiemes were all found to give similar and reasonable results, but tfie new scfiemes HLL3R and HLL5R 
were altogetfier about eigfit times more efficient in tfiis test. Tfie Roe solver-based scfieme in tfie FLASH 
code was sligfitly less dissipative, but fiad to be run at a four times lower CFL number to be stable. At RMS 
sonic Macfi number 10 only tfie new scfiemes yielded pfiysical results. We found reasonable dependence of 
dissipation on numerical resolution at tfiis Macfi number, and were able to infer a small inertial range from 
tfie velocity power spectra at 512^ resolution. Finally, we studied tfie turbulent dynamo action at RMS 
Macfi number 5 witfi tfie new scfieme. So far, tfiere fiave been very few studies of turbulent dynamos in 
tfie supersonic regime. We found dynamo- generated exponential growtfi rates of tfie magnetic pressure tfiat 
differed according to tfie type of forcing mecfianism, i.e., solenoidal versus compressive forcing [T7 1 [TS l [T9]. 
Turbulence simulations witfi tfie new scfieme across a wide range of modest to large sonic and Alfvenic Macfi 
numbers fiave been presented in [lOl [9] . 

Tfie two relaxation-based approximate Riemann solvers HLL3R and HLL5R fiave previously not been 
compared at fiigfi order and in fiigfier spatial dimensions tfian ID. In many cases we found tfiem to give 
similar results, witfi HLL3R being sligfitly more efficient. However, we presented one case, a 2D Kelvin- 
Helmfioltz instability, wfiere tfie more detailed HLL5R was significantly less viscous. Tfiis was because of 
velocity sfiears parallel to botfi tfie grid and tfie magnetic field lines. In a tfiree-dimensional turbulence run 
at Macfi 2, we also found HLL5R to be somewfiat less dissipative tfian HLL3R, giving results tfiat were very 
close to tfie Roe solver-based FLASH version. 

Tfie new FLASH MHD module is freely available upon contact witfi tfie corresponding autfior. 
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